Spatial Econometrics

Week 13

Alex Cardazzi

Old Dominion University

Topics

  • Geographic Information System Tools
    • sf, spdep
  • Spatial Econometrics
    • Moran’s I
    • Spatial Lag, Spatial Error

GIS Tools

The sf package (cheatsheet) is an extremely helpful tool for most things GIS in R.

  • Calculate distances, lengths, nearest neighbor
  • Conclude if point is inside some boundary
  • Create buffer zones, grids

The main GIS filetypes are .shp and .geojson. Be sure to know your coordinate reference system (CRS).

GIS Tools

Code 
# install.packages("sf"); install.packages("spdep")
library("sf")

# https://data.cityofnewyork.us/City-Government/Borough-Boundaries/tqmj-j8zm
nyc <- read_sf("https://data.cityofnewyork.us/api/geospatial/tqmj-j8zm?method=export&format=GeoJSON")
head(nyc)
Simple feature collection with 5 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.25559 ymin: 40.49613 xmax: -73.70001 ymax: 40.91553
Geodetic CRS:  WGS 84
# A tibble: 5 × 5
  boro_code boro_name     shape_area    shape_leng                      geometry
  <chr>     <chr>         <chr>         <chr>                 <MULTIPOLYGON [°]>
1 5         Staten Island 1623631283.36 325924.002076 (((-74.05051 40.56642, -7…
2 2         Bronx         1187189499.3  463277.240478 (((-73.89681 40.79581, -7…
3 1         Manhattan     636605816.437 359103.151368 (((-74.01093 40.68449, -7…
4 3         Brooklyn      1934169228.83 728478.125489 (((-73.86327 40.58388, -7…
5 4         Queens        3041397430.33 888238.562635 (((-73.82645 40.59053, -7…

GIS Tools

Code 
OLD_MAR <- par()$mar
par(mar = rep(0, 4))

plot(nyc$geometry)

GIS Tools

Code 
par(mar = rep(0, 4))
plot(nyc$geometry, col = alpha("tomato", .2))

GIS Tools

Code 
par(mar = rep(0, 4))
plot(nyc$geometry, border = NA, col = alpha("tomato", .2))

GIS Tools

Code 
# https://data.cityofnewyork.us/Transportation/Subway-Entrances/drex-xx56
subway <- read.csv("https://data.cityofnewyork.us/api/views/he7q-3hwy/rows.csv?accessType=DOWNLOAD")
ptz <- unlist(strsplit(gsub("POINT \\(|\\)$", "", subway$the_geom), "\\s+"))
subway$lat <- as.numeric(ptz[c(F,T)])
subway$long <- as.numeric(ptz[c(T,F)])
subway$the_geom <- NULL
head(subway)
  OBJECTID                               URL
1     1734 http://web.mta.info/nyct/service/
2     1735 http://web.mta.info/nyct/service/
3     1736 http://web.mta.info/nyct/service/
4     1737 http://web.mta.info/nyct/service/
5     1738 http://web.mta.info/nyct/service/
6     1739 http://web.mta.info/nyct/service/
                                     NAME LINE      lat      long
1 Birchall Ave & Sagamore St at NW corner  2-5 40.84917 -73.86836
2 Birchall Ave & Sagamore St at NE corner  2-5 40.84913 -73.86821
3 Morris Park Ave & 180th St at NW corner  2-5 40.84122 -73.87350
4 Morris Park Ave & 180th St at NW corner  2-5 40.84145 -73.87289
5       Boston Rd & 178th St at SW corner  2-5 40.84082 -73.87962
6  Boston Rd & E Tremont Ave at NW corner  2-5 40.84043 -73.88001

GIS Tools

Code 
subway_sf <- st_as_sf(subway, coords = c("long","lat"))
st_crs(subway_sf) <- st_crs(nyc)

GIS Tools

Code 
par(mar = rep(0, 4))
plot(nyc$geometry, border = NA, col = alpha("tomato", .2))
plot(subway_sf$geometry, pch = 20, col = alpha("black", .2), add = TRUE)

GIS Tools

Code 
head(subway_sf <- st_join(subway_sf, nyc[,"boro_name"]))
Simple feature collection with 6 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -73.88001 ymin: 40.84043 xmax: -73.86821 ymax: 40.84917
Geodetic CRS:  WGS 84
  OBJECTID                               URL
1     1734 http://web.mta.info/nyct/service/
2     1735 http://web.mta.info/nyct/service/
3     1736 http://web.mta.info/nyct/service/
4     1737 http://web.mta.info/nyct/service/
5     1738 http://web.mta.info/nyct/service/
6     1739 http://web.mta.info/nyct/service/
                                     NAME LINE boro_name
1 Birchall Ave & Sagamore St at NW corner  2-5     Bronx
2 Birchall Ave & Sagamore St at NE corner  2-5     Bronx
3 Morris Park Ave & 180th St at NW corner  2-5     Bronx
4 Morris Park Ave & 180th St at NW corner  2-5     Bronx
5       Boston Rd & 178th St at SW corner  2-5     Bronx
6  Boston Rd & E Tremont Ave at NW corner  2-5     Bronx
                    geometry
1 POINT (-73.86836 40.84917)
2 POINT (-73.86821 40.84913)
3  POINT (-73.8735 40.84122)
4 POINT (-73.87289 40.84145)
5 POINT (-73.87962 40.84082)
6 POINT (-73.88001 40.84043)

GIS Tools

Code 
par(mar = rep(0, 4))
plot(nyc$geometry, border = NA, col = alpha("tomato", .2))
plot(subway_sf$geometry, pch = 20, col = alpha("black", .2), add = TRUE)
subway_sf_buf <- st_buffer(subway_sf[1:100,], 400)
plot(subway_sf_buf$geometry, add = TRUE, col = alpha("dodgerblue", .2))

GIS Tools

Code 
st_intersects(subway_sf_buf, subway_sf)
Sparse geometry binary predicate list of length 100, where the
predicate was `intersects'
first 10 elements:
 1: 1, 2, 1878
 2: 1, 2, 1878
 3: 3, 4
 4: 3, 4
 5: 5, 6, 7, 8, 9
 6: 5, 6, 7, 8, 9
 7: 5, 6, 7, 8, 9
 8: 5, 6, 7, 8, 9
 9: 5, 6, 7, 8, 9
 10: 10, 11, 1448, 1449

GIS Tools

Code 
nn <- st_nearest_feature(subway_sf[1:7,])
# for not itself:
# st_nearest_feature(subway_sf[1:7], other_sf)
nn
[1] 2 1 4 3 6 7 6

GIS Tools

Code 
st_distance(subway_sf[1:7,], subway_sf[nn,])
Units: [m]
          [,1]      [,2]      [,3]      [,4]       [,5]       [,6]       [,7]
[1,]   12.8630    0.0000 938.99677 983.78121 1379.70011 1375.78636 1379.70011
[2,]    0.0000   12.8630 939.80671 985.05959 1385.08336 1381.02949 1385.08336
[3,]  985.0596  983.7812  57.10861   0.00000  554.28722  541.52211  554.28722
[4,]  939.8067  938.9968   0.00000  57.10861  608.99494  596.54355  608.99494
[5,] 1332.5241 1327.0734 570.65182 517.15751   53.17395   54.21955   53.17395
[6,] 1385.0834 1379.7001 608.99494 554.28722    0.00000   16.98492    0.00000
[7,] 1381.0295 1375.7864 596.54355 541.52211   16.98492    0.00000   16.98492

GIS Tools

Code 
st_distance(subway_sf[1:7,], subway_sf[nn,], by_element = TRUE)
Units: [m]
[1] 12.86300 12.86300 57.10861 57.10861 53.17395 16.98492 16.98492

Tobler’s First Law of Geography

“Everything is related to everything else, but near things are more related than distant things.”

We often need to consider spatial interactions and spillovers in our models. Sometimes because they are interesting, other times because we’d otherwise violate assumptions.

  • Spread of Disease
  • Effects of Pollution
  • Effects of legal marijuana

Modifiable Areal Unit Problem (MAUP)

MAUP is a source of statistical bias that can result when point-based measures are aggregated into larger areas.

  • Effect of luminosity on crime:

    - Crime Counts? Choose an area (corner, block, neighborhood (?), county, state)
- Whether a crime happens? Need all possible criminal interactions
- Use CI (and counterfactuals) to inform this

  • Effect of pavement smoothness on crash rates

    - What units should we use for this?
- A mile? A kilometer? Distance between exits? Routes? Route-by-County? A grid?

Modifiable Areal Unit Problem (MAUP)

Autocorrelation

Classically, autocorrelation is defined as when nearby observations are correlated with one another.

  • \(\text{cor}(Y_t, Y_{t-1})\)

We can generalize the word “nearby” to mean temporally and spatially.

  • \(\text{cor}(Y_{i,t}, Y_{i, t-1})\)
  • \(\text{cor}(Y_{i,t}, Y_{j, t})\)

Autocorrelation

Now, we can think of cross sectional (and panel) data the way we thought of time series data.

  • OLS: \(Y_i = f(X_i, e_i)\)
  • Spatial lag: \(Y_i = f(X_i, X_j, e_i)\) or \(Y_i = f(X_i, Y_j, e_i)\)
  • Spatial error: \(Y_i = f(X_i, e_i, e_j)\)
  • Spatial heterogeneity: \(Y_i = f(X_i, e_i, \gamma)\)

Spatial Autocorrelation

Spatial Weights

To calculate spatial autocorrelation, define “nearby”:

  • Continuous Distance (m meters away)
  • Binary Distance (within some m meters)
  • Continguity (next to, touching)
  • Interactions
  • K-nearest

Your results will be strongly dependent on how you define “nearby”.

Spatial Weights

For simplicity, we will only consider rook and queen contiguity.

Spatial Weights

Code 
# https://data.cityofnewyork.us/City-Government/2010-Neighborhood-Tabulation-Areas-NTAs-/cpf4-rkhq
nta <- read_sf("https://data.cityofnewyork.us/api/geospatial/cpf4-rkhq?method=export&format=GeoJSON")
head(nta)
Simple feature collection with 6 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -73.98843 ymin: 40.6087 xmax: -73.74359 ymax: 40.88596
Geodetic CRS:  WGS 84
# A tibble: 6 × 8
  ntacode shape_area    county_fips ntaname              shape…¹ boro_…² boro_…³
  <chr>   <chr>         <chr>       <chr>                <chr>   <chr>   <chr>  
1 QN08    77412747.7604 081         St. Albans           45401.… Queens  4      
2 BK69    20528197.0577 047         Clinton Hill         23971.… Brookl… 3      
3 BK46    17782095.6686 047         Ocean Parkway South  21975.… Brookl… 3      
4 BX28    25666124.7136 005         Van Cortlandt Villa… 21945.… Bronx   2      
5 QN55    82461393.5232 081         South Ozone Park     36708.… Queens  4      
6 BK40    14041667.948  047         Windsor Terrace      19033.… Brookl… 3      
# … with 1 more variable: geometry <MULTIPOLYGON [°]>, and abbreviated variable
#   names ¹​shape_leng, ²​boro_name, ³​boro_code

Spatial Weights

Code 
par(mar = rep(0, 4))
plot(nta$geometry, col = alpha(nta$boro_code, .6), border = "white")

Spatial Weights

Code 
par(mar = rep(0, 4))
nta$shape_area <- as.numeric(nta$shape_area)
colz <- (nta$shape_area - min(nta$shape_area)) / (max(nta$shape_area) - min(nta$shape_area))
plot(nta$geometry, col = alpha(nta$boro_code, colz), border = "white")

Spatial Weights

Code 
par(mar = rep(0, 4))
plot(nta$geometry, col = alpha(nta$boro_code, .6), border = "white")
plot(subway_sf$geometry, col = alpha("black", .2), pch = 20, add = TRUE)

Spatial Weights

Code 
library("spdep")

nb_r <- poly2nb(nta, queen = FALSE)
nb_q <- poly2nb(nta, queen = TRUE)

Spatial Weights

Code 
summary(nb_r)
Neighbour list object:
Number of regions: 195 
Number of nonzero links: 1046 
Percentage nonzero weights: 2.750822 
Average number of links: 5.364103 
1 region with no links:
45
Link number distribution:

 0  1  2  3  4  5  6  7  8  9 10 17 24 27 37 
 1  1  9 22 41 56 33 16  8  2  2  1  1  1  1 
1 least connected region:
82 with 1 link
1 most connected region:
42 with 37 links
Code 
summary(nb_q)
Neighbour list object:
Number of regions: 195 
Number of nonzero links: 1094 
Percentage nonzero weights: 2.877055 
Average number of links: 5.610256 
1 region with no links:
45
Link number distribution:

 0  1  2  3  4  5  6  7  8  9 10 17 24 27 38 
 1  1  8 19 35 55 33 20 12  5  2  1  1  1  1 
1 least connected region:
82 with 1 link
1 most connected region:
42 with 38 links

Spatial Weights

Code 
par(mar = rep(0, 4))
coords <- st_centroid(nta$geometry)
plot(nta$geometry, col = alpha(nta$boro_code, .6), border = "white")
plot(nb_r, coords, pch = 20, add = TRUE, col = alpha("tomato", .6))

Spatial Weights

Code 
par(mar = rep(0, 4))
plot(nta$geometry, col = alpha(nta$boro_code, .6), border = "white")
plot(nb_q, coords, pch = 20, add = TRUE, col = alpha("tomato", .6))

Spatial Weights

Code 
set.ZeroPolicyOption(TRUE) # get.ZeroPolicyOption()
w <- nb2listw(nb_q, style = "W")
Code 
w$neighbours[1:4]
[[1]]
[1]  26  28  29  42  46  59 120 141 193

[[2]]
[1]  15  24  36  54  61 138

[[3]]
[1]  35  40  70  96 145

[[4]]
[1] 137 167 170 171 179
Code 
w$weights[1:4]
[[1]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[2]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[3]]
[1] 0.2 0.2 0.2 0.2 0.2

[[4]]
[1] 0.2 0.2 0.2 0.2 0.2

Spatial Autocorrelation

Code 
hous <- read_sf("https://data.cityofnewyork.us/api/geospatial/kyz5-72x5?method=export&format=GeoJSON")
m <- match(nta$ntacode, hous$nta2010)
nta$units <- as.numeric(hous$cenunits10[m])
nta$permits <- as.numeric(hous$filed[m])
nta$permits <- ifelse(is.na(nta$permits), 0, nta$permits)

tmp <- st_join(subway_sf, nta[,"ntacode"])
tmp <- as.data.frame(table(tmp$ntacode))
nta$subway <- tmp$Freq[match(nta$ntacode, tmp$Var1)]
nta$subway <- ifelse(is.na(nta$subway), 0, nta$subway)

head(nta[,c(4, 6, 9:11)])
Simple feature collection with 6 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -73.98843 ymin: 40.6087 xmax: -73.74359 ymax: 40.88596
Geodetic CRS:  WGS 84
# A tibble: 6 × 6
  ntaname               boro_name units permits subway                  geometry
  <chr>                 <chr>     <dbl>   <dbl>  <dbl>        <MULTIPOLYGON [°]>
1 St. Albans            Queens    16170      36      0 (((-73.75205 40.70523, -…
2 Clinton Hill          Brooklyn  16997     417     17 (((-73.95337 40.68064, -…
3 Ocean Parkway South   Brooklyn   8065      -6      5 (((-73.97075 40.62563, -…
4 Van Cortlandt Village Bronx     18670     343      6 (((-73.88705 40.88435, -…
5 South Ozone Park      Queens    23011      36     11 (((-73.80577 40.68293, -…
6 Windsor Terrace       Brooklyn   9684      17      7 (((-73.98017 40.66115, -…

Moran’s I

We want to come up with a quantity similar to a (pearson) correlation coefficient that runs from -1 to 1.

Let’s use correlation as a motivation:

\[r = \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum {(x_i - \bar{x})^2}\sum {(y_i - \bar{y})^2}}}\]

\[I = \frac{N \sum_i \sum_j w_{i,j} (x_i - \bar{x})(x_j - \bar{x})}{\sum_i \sum_j w_{i,j} \sum_i (x_i - \bar{x})^2} \]

Moran’s I

Moran’s I

Using a Queen contiguity weight matrix:

\[\bar{x} = \frac{1}{2}\]

\[\sum w_{i,j} = 32\]

\[\sum (x_i - \bar{x})^2 = (1-\frac{1}{2})^2 + ... + (0-\frac{1}{2})^2 = 8 * \frac{1}{4} = 2\]

\[\sum_i \sum_j w_{i,j} (x_i - \bar{x})(x_j - \bar{x}) = 3 (1 * \frac{1}{2} * \frac{1}{2}) + 5 (0 * \frac{1}{2} * \frac{1}{2}) + ... = 4\]

\[I_{\text{top}} = \frac{8 * 4}{32 * 2} = \frac{1}{2}\] \[I_{\text{bottom}} = \frac{8 * -2}{32 * 2} = -\frac{1}{4}\]

Moran’s I

Code 
moran(nta$permits, w, n = nrow(nta), S0 = Szero(w))
moran.test(nta$permits, w)
$I
[1] 0.1081634

$K
[1] 26.6402


    Moran I test under randomisation

data:  nta$permits  
weights: w  n reduced by no-neighbour observations
  

Moran I statistic standard deviate = 2.712, p-value = 0.003344
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.107608672      -0.005181347       0.001729663

Moran’s I

Code 
moran.plot(nta$permits, w)

Spatial Weights

Code 
perms.lag <- lag.listw(w, nta$permits)
subwy.lag <- lag.listw(w, nta$subway)

Spatial Weights

Code 
par(mfrow = c(1, 2))
plot(nta$permits, perms.lag)
plot(nta$subway, subwy.lag)

Modeling Change in Building Permits

Suppose we want to model the change in building permits as a function of public transportation accessability. If there is more easily accessible transportation, perhaps demand (housing stock) will grow faster. Our model would look like the following:

\[\text{Permits}_i = \alpha + \beta*\text{Subway}_i + e_i\]

Modeling Change in Building Permits

Code 
par(mfrow = c(1, 2), mar = rep(0, 4))
plot(nta$geometry, col = alpha("tomato", asinh(nta$subway)/max(asinh(nta$subway))))
plot(nta$geometry, col = alpha("tomato", asinh(nta$permits)/max(asinh(nta$permits))))

Modeling Change in Building Permits

Code 
r <- lm(permits ~ subway, nta)
summary(r)

Call:
lm(formula = permits ~ subway, data = nta)

Residuals:
   Min     1Q Median     3Q    Max 
-874.8 -170.9 -129.8  115.7 2680.2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  174.796     28.666   6.098 5.75e-09 ***
subway         8.098      1.534   5.279 3.49e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 339.7 on 193 degrees of freedom
Multiple R-squared:  0.1262,    Adjusted R-squared:  0.1216 
F-statistic: 27.86 on 1 and 193 DF,  p-value: 3.485e-07

Modeling Change in Building Permits

Code 
lm.morantest(r, w)

    Global Moran I for regression residuals

data:  
model: lm(formula = permits ~ subway, data = nta)
weights: w

Moran I statistic standard deviate = 2.4555, p-value = 0.007035
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.101818703     -0.006759151      0.001955277

Modeling Change in Building Permits

Code 
color_func <- colorRampPalette(c("dodgerblue", "tomato"))
colorz <- color_func(length(unique(r$residuals)))
par(mar = rep(0, 4))
plot(nta$geometry, col = alpha(colorz[order(r$residuals)], .6))

Spatial Lag in X

Good news!

A model with spatially lagged x variables does not impact the OLS assumptions! You can include these in the model without any complications.

Imagine a model of crime as a function of unemployment. You could model this as the following:

\[\text{Crime}_i = \alpha + \beta_1 \text{Unemployment}_i + \beta_2 \text{Unemployment}_j + e_i\]

Spatial Lag in X

Code 
r1 <- spatialreg::lmSLX(permits ~ subway,
                        data = nta, w,
                        zero.policy = TRUE)
summary(r1)

Call:
lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
    data = as.data.frame(x), weights = weights)

Residuals:
   Min     1Q Median     3Q    Max 
-893.0 -172.4 -127.8  114.0 2670.8 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 180.2075    34.1184   5.282 3.45e-07 ***
subway        8.3230     1.7169   4.848 2.57e-06 ***
lag.subway   -0.7365     2.5040  -0.294    0.769    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 340.5 on 192 degrees of freedom
Multiple R-squared:  0.1266,    Adjusted R-squared:  0.1175 
F-statistic: 13.91 on 2 and 192 DF,  p-value: 2.285e-06

Spatial Error

If you are still getting spatially correlated residuals, there might be some unobserved, spatially correlated factor in the error term. Imagine modeling average income in a county as a function of average education. While location/distance might not influence income or education, the way counties are created might make closer counties more similar to nearby counties than far away counties. This suggests that the error in the model is correlated across space.

Importantly, a spatial error model only induces inefficient parameter estimates (bad standard errors).

Spatial Error

\[\begin{align} \text{Income}_i &= X_i\beta + \epsilon_i \\ \epsilon_i &= \lambda w_i\epsilon_j + e_i\end{align}\]

Code 
spatialreg::errorsarlm(permits ~ subway,
                        data = nta, w,
                        zero.policy = TRUE) -> r2
summary(r2)

Call:spatialreg::errorsarlm(formula = permits ~ subway, data = nta, 
    listw = w, zero.policy = TRUE)

Residuals:
    Min      1Q  Median      3Q     Max 
-886.16 -169.25 -120.44  119.76 2698.23 

Type: error 
Regions with no neighbours included:
 45 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) 170.3031    33.9190  5.0209 5.144e-07
subway        8.3914     1.5910  5.2744 1.332e-07

Lambda: 0.21156, LR test value: 4.2897, p-value: 0.038344
Asymptotic standard error: 0.10441
    z-value: 2.0262, p-value: 0.042741
Wald statistic: 4.1056, p-value: 0.042741

Log likelihood: -1409.994 for error model
ML residual variance (sigma squared): 110760, (sigma: 332.81)
Number of observations: 195 
Number of parameters estimated: 4 
AIC: 2828, (AIC for lm: 2830.3)

Spatial Lag

Imagine a model of median income by county. If earnings increase in Norfolk, those residents might go spend more money in Virginia Beach, thus increasing earnings for Virginia Beach. This would be positive spatial autocorrelation, or spillover.

Imagine a model of unemployment by county. If unemployment decreases in Norfolk (and total jobs remains the same) it could be that they are cannibalizing jobs from Virginia Beach residents. This would be negative spillover.

Spatial Lag

\[\begin{align} \text{Income}_i &= X_i\beta + \epsilon_i \\ \epsilon_i &= \rho w_i y_j + e_i\end{align}\]

Code 
spatialreg::lagsarlm(permits ~ subway,
                        data = nta, w,
                        zero.policy = TRUE) -> r3
summary(r3)

Call:spatialreg::lagsarlm(formula = permits ~ subway, data = nta, 
    listw = w, zero.policy = TRUE)

Residuals:
    Min      1Q  Median      3Q     Max 
-838.20 -160.26 -120.47  113.43 2717.32 

Type: lag 
Regions with no neighbours included:
 45 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) 130.2620    35.9966  3.6187 0.0002961
subway        7.8174     1.5383  5.0819 3.736e-07

Rho: 0.18002, LR test value: 3.5057, p-value: 0.061158
Asymptotic standard error: 0.099035
    z-value: 1.8178, p-value: 0.069098
Wald statistic: 3.3043, p-value: 0.069098

Log likelihood: -1410.386 for lag model
ML residual variance (sigma squared): 111470, (sigma: 333.88)
Number of observations: 195 
Number of parameters estimated: 4 
AIC: 2828.8, (AIC for lm: 2830.3)
LM test for residual autocorrelation
test value: 0.35178, p-value: 0.55311
Code 
spatialreg::impacts(r3, listw = w)
Impact measures (lag, exact):
         Direct Indirect    Total
subway 7.866091 1.667634 9.533724

Final Project

Final Project

Next week: “brown bag” presentations

  • You will be graded on your presentation but it does not need to be perfect. You will get feedback from myself as well as at least one other person. Please implement this feedback into your final submission
  • Your submission should be a memo + manual (like before).
    • Why is this important?
    • Why/how did you select the variables you did?
    • Why/how did you select the functional form?
    • How does your model preform?
    • What is your forecast?
    • Do you trust it?

Two weeks: projects are be due Dec. 14 at 12:01 am.

ECON 707/807: Econometrics II